module studentscr
implicit none
contains
subroutine stllik(m,theta,likresm)
use linear_operators
use matrixfuns
use datahadler
integer, intent(in):: m
double precision, intent(in):: theta(m)
double precision, intent(out):: likresm
integer, parameter:: a=1000,b=3
integer T,n,i,j
double precision y(a,b)
double precision, parameter:: c9=.999d0
common T,n,y
double precision, parameter:: pi=3.14159265358979d0
double precision lamb(T),delta(n),phi(n),c(n),phi0(n),phi1(n),phi2(n),&
                &alpha1,alpha2,eta,llik,sigmat(n,n),gamma(n),gammat(n,n),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n),alpha0&
				&,fttm,wttm,jit,cgt,cfun
!call exportvec(theta,"thstint.txt",5)


delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall
eta=theta(2*n+6)

forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


gamma=phi0/(1.0d0-phi1-phi2)


gammat=diag(gamma)

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)

igammat=diag(1.0d0/gamma)
sigmat=(cc*lamb(1))+gammat
isigmat=igammat-lamb(1)*(((igammat.x.cc).x.igammat)&
     &/(lamb(1)*dot_product(c,matmul(igammat,c))+1.0d0))



jit=dot_product(epst(1,:),matmul(isigmat,epst(1,:)))
call cgive(eta)
cgt=coefg(jit)
llik=cfun-.5d0*dlog(det(sigmat))+cgt


do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0&
	       &+alpha1*(((fttm)**2.0d0)+wttm)+alpha2*lamb(i-1)


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)

	igammat=diag(1.0d0/gamma)

	sigmat=(cc*lamb(i))+gammat
	isigmat=igammat-lamb(i)*(((igammat.x.cc).x.igammat)&
	     &/(lamb(i)*dot_product(c,matmul(igammat,c))+1.0d0))

	jit=dot_product(epst(i,:),matmul(isigmat,epst(i,:)))
	
	cgt=coefg(jit)
	llik=llik+cfun-.5d0*dlog(det(sigmat))+cgt

end do
likresm=-llik
contains
subroutine cgive(eta)
double precision g1,g2,s1,s2,eta
external dlgams
if (eta<8.0d-4) then
	cfun=-.5d0*n*dlog(2*pi)+.25d0*n*(n+2.0d0)*eta-(n*(n+2)*(n-5)/13.0d0)*eta*eta&
	&+(n*(n+2)*(n*n-6*n+16)/24.0d0)*eta*eta*eta
else
	call dlgams((n*eta+1.0d0)/(2.0d0*eta),g1,s1)
	call dlgams(1.0d0/(2.0d0*eta),g2,s2)
	cfun=g1-g2-.5d0*n*dlog((1.0d0-2.0d0*eta)/eta)-.5d0*n*dlog(pi)
end if
end subroutine cgive

function coefg(ji)
double precision ji,coefg
    if (eta>.03d0) then
	    coefg=-((n*eta+1.0d0)/(2.0d0*eta))*dlog(1.0d0+(eta/(1.0d0-2.0d0*eta))*ji)
    elseif ((ji*eta)>.001d0) then
        coefg=-((n*eta+1.0d0)/(2.0d0*eta))*dlog(1.0d0+(eta/(1.0d0-2.0d0*eta))*ji)
    else
        coefg=-.5d0*ji+(-.5d0*(n+2d0)*ji+.25*ji**2.0d0)*eta&
		&+.5*(-2.0d0*(n+2.0d0)*ji+.5d0*(n+4.0d0)*ji**2.0d0-(1.0d0/3.0d0)*ji**3.0d0)*eta**2.0d0&
           &+(1.0d0/6.0d0)*(-13.0d0*(n+2.0d0)*ji+6.0d0*(n+3.0d0)*ji**2.0d0-(n+6.0d0)*ji**3.0d0+(1.0d0/8.0d0)*ji**4)*eta**3.0d0 &
           &+(1.0/24.0d0)*(-96.0d0*(n+2.0d0)*ji+24.0d0*(8.0d0+3.0d0*n)*ji**2.0d0-24.0d0*(n+4.0d0)*ji**3.0d0+3.0d0*(n+8.0d0)*ji**4.0d0-(13.0d0/5.0d0)*ji**5.0d0)*eta**4.0d0 &
           &+(1.0d0/130.0d0)*(-960.0d0*(n+2.0d0)*ji+600.0d0*(2.0d0*n+5.0d0)*ji**2.0d0-1440.0d0*(3.0d0*n+10.0d0)*ji**3.0d0+130.0d0*(n+5.0d0)*ji**4.0d0-13.0d0*(n+10.0d0)*ji**5.0d0+10.0d0*ji**6.0d0)*eta**5.0d0
    end if
end function coefg
end subroutine stllik


subroutine stscr(m,theta,grad)
!lambda=1
use linear_operators
use matrixfuns
use datahadler
integer, intent(in):: m
double precision, intent(in):: theta(m)
double precision, intent(out):: grad(m)
integer, parameter:: a=1000,b=3
integer T,n,i,j,ii
double precision y(a,b)
double precision, parameter:: c9=.999d0
common T,n,y
double precision lamb(T),delta(n),c(n),phi0(n),phi1(n),phi2(n),eta,gamma(n),&
                &sigmat(n,n),gammat(n,n),alpha1,alpha2,&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n)
double precision ddelta(m,n),dc(m,n),dgamma(m,n),ent(n,n*n)&
               &,dmut(m,n),dalpha1(m),dalpha2(m),alpha0,dalpha0(m),deta(m)&
			   &,dlambda(m),ident(n,n),sigeet(n),siget(n),wttm,fttm&
			   &,dfttm(m),dwttm(m),matwtt(m),metn(n),dphi0(m,n),dphi1(m,n),dphi2(m,n)&
			   &,jit,dcfun
		
ent=matent(n)
ident=eye(n)
delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall
eta=theta(2*n+6)


forall(i=1:m,j=1:n)
	ddelta(i,j)=0.0d0
	dc(i,j)=0.0d0
	dgamma(i,j)=0.0d0
	dphi0(i,j)=0.0d0
	dphi1(i,j)=0.0d0
	dphi2(i,j)=0.0d0
	dmut(i,j)=0.0d0
end forall
forall(i=1:n)
	ddelta(i,i)=1.0d0
	dphi0(2*n+1,i)=1.0d0
	dphi1(2*n+4,i)=c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4))
	dphi2(2*n+4,i)=-(c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4)))*(dsin(theta(2*n+5))**2.0d0)
	dphi2(2*n+5,i)=(c9-phi1(i))*2.0d0*dsin(theta(2*n+5))*dcos(theta(2*n+5))
end forall
forall(i=2:n)
	dc(n+i-1,i)=1.0d0
end forall

forall(i=1:m)
	dalpha0(i)=0.0d0
	dalpha1(i)=0.0d0
	dlambda(i)=0.0d0
	deta(i)=0.0d0
end forall
dalpha2=dalpha1
dalpha0(2*n)=1.0d0
dalpha1(2*n+2)=c9*2.0d0*dsin(theta(2*n+2))*dcos(theta(2*n+2))
dalpha2(2*n+2)=-dalpha1(2*n+2)*(dsin(theta(2*n+3))**2.0d0)
dalpha2(2*n+3)=(c9-alpha1)*2.0d0*dsin(theta(2*n+3))*dcos(theta(2*n+3))
deta(2*n+6)=1.0d0

forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:n)
	dmut(i,i)=1.0d0
end forall

gamma=phi0/(1.0d0-phi1-phi2)

gammat=diag(gamma)
igammat=diag(1.0d0/gamma)
forall (i=1:m,j=1:n)
	dgamma(i,j)=((1.0d0-phi1(j)-phi2(j))*dphi0(i,j)-phi0(j)*(-dphi1(i,j)-dphi2(i,j)))&
	           &/(1.0d0-phi1(j)-phi2(j))**2.0d0
end forall
dlambda=((1.0d0-alpha1-alpha2)*dalpha0-alpha0*(-dalpha1-dalpha2))/(1.0d0-alpha1-alpha2)**2.0d0

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)
sigmat=lamb(1)*cc+gammat
isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(1))))

call dcgive(eta)
jit=dot_product(epst(1,:),matmul(isigmat,epst(1,:)))

forall(i=1:n)
	sigeet(i)=((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*(dot_product(isigmat(i,:),epst(1,:))*&
	          & dot_product(isigmat(:,i),epst(1,:)))-isigmat(i,i)
end forall
siget=matmul(isigmat,epst(1,:))

grad=((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*matmul(dmut,siget)&
    &+.5d0*matmul(dgamma,sigeet)&
	&+.5d0*dlambda*(((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*(dot_product(c,siget)**2.0)-dot_product(c,matmul(isigmat,c)))&
	&+lamb(1)*(((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*dot_product(c,siget)*matmul(dc,siget)&
	&-matmul(matmul(dc,isigmat),c))+(dcfun+dcoefg(jit))*deta



!t>1
do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0+alpha1*((fttm**2.0d0)+wttm)+alpha2*lamb(i-1)

	forall(ii=1:n)
	    metn(ii)=dot_product(igammat(ii,:),c)*dot_product(igammat(:,ii),epst(i-1,:))
	end forall
	matwtt=matmul(matmul(dgamma,ent),vec2((igammat.x.cc).x.igammat))

	dwttm=-(wttm*wttm)*(-(lamb(i-1)**(-2.0))*dlambda+2.0d0*matmul(dc,matmul(igammat,c))&
	      &-matwtt)
	
	dfttm=dwttm*dot_product(c,matmul(igammat,epst(i-1,:)))&
	     &+wttm*matmul(dc,matmul(igammat,epst(i-1,:)))&
		 &-wttm*matmul(dgamma,metn)-wttm*matmul(dmut,matmul(igammat,c))


	forall (ii=1:m,j=1:n)
		dgamma(ii,j)=dphi0(ii,j)+dphi1(ii,j)*(((epst(i-1,j)-c(j)*fttm)**2.0d0)&
		                                                                   &+(c(j)*c(j))*wttm)&
		           &+dphi2(ii,j)*gamma(j)&
				   &+phi1(j)*2.0d0*(epst(i-1,j)-c(j)*fttm)*&
				   &(-dmut(ii,j)-dc(ii,j)*fttm-c(j)*dfttm(ii))&
				   &+phi1(j)*2.0d0*c(j)*dc(ii,j)*wttm+phi1(j)*c(j)*c(j)*dwttm(ii)+phi2(j)*dgamma(ii,j)
	end forall


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)
	igammat=diag(1.0d0/gamma)


	sigmat=lamb(i)*cc+gammat
	isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(i))))

	jit=dot_product(epst(i,:),matmul(isigmat,epst(i,:)))

	forall(ii=1:n)
		sigeet(ii)=((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*(dot_product(isigmat(ii,:),epst(i,:))*&
		          &   dot_product(isigmat(:,ii),epst(i,:)))-isigmat(ii,ii)
	end forall
	siget=matmul(isigmat,epst(i,:))

	dlambda=dalpha0+dalpha1*((fttm*fttm)+wttm)+dalpha2*(lamb(i-1))&
	      &+alpha1*(2.0d0*fttm*dfttm+dwttm)+alpha2*dlambda

	grad=grad+((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*matmul(dmut,siget)&
	    &+.5d0*matmul(dgamma,sigeet)&
		&+.5d0*dlambda*(((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*(dot_product(c,siget)**2.0)-dot_product(c,matmul(isigmat,c)))&
		&+lamb(i)*(((N*eta+1.0d0)/(1.0d0-2.0d0*eta+eta*jit))*dot_product(c,siget)*matmul(dc,siget)&
		&-matmul(matmul(dc,isigmat),c))+(dcfun+dcoefg(jit))*deta

end do

grad=-grad
contains
subroutine dcgive(eta)
double precision eta,dpsi

if (eta<8.0d-4) then
	dcfun=.25d0*n*(n+2.0d0)-2.0d0*(n*(n+2)*(n-5)/13.0d0)*eta&
	     &+3.0d0*(n*(n+2)*(n*n-6*n+16)/24.0d0)*eta*eta
else
    dcfun=(n/(2.0d0*eta*(1.0d0-2.0d0*eta)))&
	    &-(1.0d0/(2.0d0*eta**2.0d0))*(dpsi((n*eta+1.0d0)/(2.0d0*eta))-dpsi(1.0d0/(2.0d0*eta)))
end if
end subroutine dcgive

function dcoefg(ji)
double precision ji,dcoefg
    if (eta>.03d0) then
        dcoefg=-((n*eta+1.0d0)/(2.0d0*eta*(1.0d0-2.0d0*eta)))*(ji/(1.0d0-2.0d0*eta+eta*ji))&
		&+(1.0d0/(2.0d0*eta**2.0d0))*dlog(1.0d0+(eta/(1.0d0-2.0d0*eta))*ji)		
    elseif ((ji*eta)>.001d0) then
        dcoefg=-((n*eta+1.0d0)/(2.0d0*eta*(1.0d0-2.0d0*eta)))*(ji/(1.0d0-2.0d0*eta+eta*ji))&
		&+(1.0d0/(2.0d0*eta**2.0d0))*dlog(1.0d0+(eta/(1.0d0-2.0d0*eta))*ji)
    else
        dcoefg=(-.5d0*(n+2d0)*ji+.25*ji**2.0d0)&
		&+(-2.0d0*(n+2.0d0)*ji+.5d0*(n+4.0d0)*ji**2.0d0-(1.0d0/3.0d0)*ji**3.0d0)*eta&
           &+(1.0d0/6.0d0)*(-13.0d0*(n+2.0d0)*ji+6.0d0*(n+3.0d0)*ji**2.0d0-(n+6.0d0)*ji**3.0d0+(1.0d0/8.0d0)*ji**4)*3.0d0*eta**2.0d0 &
           &+(1.0/24.0d0)*(-96.0d0*(n+2.0d0)*ji+24.0d0*(8.0d0+3.0d0*n)*ji**2.0d0-24.0d0*(n+4.0d0)*ji**3.0d0+3.0d0*(n+8.0d0)*ji**4.0d0-(13.0d0/5.0d0)*ji**5.0d0)*4.0d0*eta**3.0d0 &
           &+(1.0d0/130.0d0)*(-960.0d0*(n+2.0d0)*ji+600.0d0*(2.0d0*n+5.0d0)*ji**2.0d0-1440.0d0*(3.0d0*n+10.0d0)*ji**3.0d0+130.0d0*(n+5.0d0)*ji**4.0d0-13.0d0*(n+10.0d0)*ji**5.0d0+10.0d0*ji**6.0d0)*5.0d0*eta**4.0d0
    end if
end function dcoefg

end subroutine stscr

subroutine stsimul(m,T,n,theta,y)
use linear_operators
use matrixfuns
use ghsrnd
integer, intent(in):: m,T,n
double precision, intent(in):: theta(m)
double precision, intent(out):: y(T,n)
integer i,j,irank
double precision, parameter:: c9=.999d0
double precision, parameter:: pi=3.14159265358979d0
double precision lamb(T),delta(n),phi(n),c(n),phi0(n),phi1(n),phi2(n),&
                &alpha1,alpha2,eta,llik,sigmat(n,n),gamma(n),gammat(n,n),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n)&
				&,fttm,wttm,rsigmat(n,n),epstar(T,n),beta(n),tol,dmach,alpha0
delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)


forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall
eta=theta(2*n+6)

forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:n)
	beta(i)=0.0d0
end forall
epstar=vsktdrnd(eta,beta,T,n)

gamma=phi0/(1.0d0-phi1-phi2)

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)

gammat=diag(gamma)
sigmat=cc*lamb(1)+gammat

tol = 100.0*dmach(4)

call dchfac(n,sigmat,n,tol,irank,rsigmat,n)
rsigmat=.t.rsigmat

epst(1,:)=matmul(rsigmat,epstar(1,:))
y(1,:)=delta+epst(1,:)


do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0+alpha1*(((fttm)**2.0d0)+wttm)+alpha2*lamb(i-1)


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)

	sigmat=(cc*lamb(i))+gammat

	call dchfac(n,sigmat,n,tol,irank,rsigmat,n)
	rsigmat=.t.rsigmat

	epst(i,:)=matmul(rsigmat,epstar(i,:))
	y(i,:)=delta+epst(i,:)
end do

end subroutine stsimul

!!!EM algorithm

subroutine max_step(m,theta,likresm)
use linear_operators
use matrixfuns
use datahadler
integer, intent(in):: m
double precision, intent(in):: theta(m)
double precision, intent(out):: likresm
integer, parameter:: a=1000,b=3
integer T,n,i,j
double precision y(a,b),eta_0
double precision, parameter:: c9=.999d0
common T,n,y
common /eta_prev_step/ eta_0
double precision, parameter:: pi=3.14159265358979d0
double precision lamb(T),delta(n),phi(n),c(n),phi0(n),phi1(n),phi2(n),&
                &alpha1,alpha2,eta,llik,sigmat(n,n),gamma(n),gammat(n,n),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n),alpha0&
				&,fttm,wttm,jit,cgt,cfun,nu_0,nu
double precision:: dpsi,e_xi,e_ixi,e_logxi,g1,s1

!!!
nu_0=1.0d0/eta_0
e_xi=nu_0
e_ixi=1.0d0/(nu_0-2.0d0)
e_logxi=dlog(2.0d0)+dpsi(.5d0*nu_0)

!!!

delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall
eta=theta(2*n+6)
nu=1.0d0/eta

forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


gamma=phi0/(1.0d0-phi1-phi2)


gammat=diag(gamma)

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)

igammat=diag(1.0d0/gamma)
sigmat=(cc*lamb(1))+gammat
isigmat=igammat-lamb(1)*(((igammat.x.cc).x.igammat)&
     &/(lamb(1)*dot_product(c,matmul(igammat,c))+1.0d0))



jit=dot_product(epst(1,:),matmul(isigmat,epst(1,:)))

call dlgams(1.0d0/(2.0d0*eta),g1,s1)

llik=-.5d0*n*dlog(2.0d0*pi*(nu-2.0d0))+.5d0*n*e_logxi&
&-.5d0*dlog(det(sigmat))&
&-.5d0*jit*(e_xi/(nu-2.0d0))&
&-(1.0d0/(2.0d0*eta))*dlog(2.0d0)-g1+((1.0d0/(2.0d0*eta))-1.0d0)*e_logxi-.5d0*e_xi


do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0&
	       &+alpha1*(((fttm)**2.0d0)+wttm)+alpha2*lamb(i-1)


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)

	igammat=diag(1.0d0/gamma)

	sigmat=(cc*lamb(i))+gammat
	isigmat=igammat-lamb(i)*(((igammat.x.cc).x.igammat)&
	     &/(lamb(i)*dot_product(c,matmul(igammat,c))+1.0d0))

	jit=dot_product(epst(i,:),matmul(isigmat,epst(i,:)))
	
	llik=llik-.5d0*n*dlog(2.0d0*pi*(nu-2.0d0))+.5d0*n*e_logxi&
	&-.5d0*dlog(det(sigmat))&
	&-.5d0*jit*(e_xi/(nu-2.0d0))&
	&-(1.0d0/(2.0d0*eta))*dlog(2.0d0)-g1+((1.0d0/(2.0d0*eta))-1.0d0)*e_logxi-.5d0*e_xi

end do
likresm=-llik
end subroutine max_step

subroutine max_step_scr(m,theta,grad)
!lambda=1
use linear_operators
use matrixfuns
use datahadler
integer, intent(in):: m
double precision, intent(in):: theta(m)
double precision, intent(out):: grad(m)
integer, parameter:: a=1000,b=3
integer T,n,i,j,ii
double precision y(a,b),eta_0
double precision, parameter:: c9=.999d0
common T,n,y
common /eta_prev_step/ eta_0
double precision lamb(T),delta(n),c(n),phi0(n),phi1(n),phi2(n),eta,gamma(n),&
                &sigmat(n,n),gammat(n,n),alpha1,alpha2,&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n)
double precision ddelta(m,n),dc(m,n),dgamma(m,n),ent(n,n*n)&
               &,dmut(m,n),dalpha1(m),dalpha2(m),alpha0,dalpha0(m),deta(m)&
			   &,dlambda(m),ident(n,n),sigeet(n),siget(n),wttm,fttm&
			   &,dfttm(m),dwttm(m),matwtt(m),metn(n),dphi0(m,n),dphi1(m,n),dphi2(m,n)&
			   &,jit,dcfun
double precision:: dpsi,e_xi,e_ixi,e_logxi,g1,s1,nu_0,dlnui,nu

!!!
nu_0=1.0d0/eta_0
e_xi=nu_0
e_ixi=1.0d0/(nu_0-2.0d0)
e_logxi=dlog(2.0d0)+dpsi(.5d0*nu_0)

!!!
		
ent=matent(n)
ident=eye(n)
delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall
eta=theta(2*n+6)
nu=1.0d0/eta


forall(i=1:m,j=1:n)
	ddelta(i,j)=0.0d0
	dc(i,j)=0.0d0
	dgamma(i,j)=0.0d0
	dphi0(i,j)=0.0d0
	dphi1(i,j)=0.0d0
	dphi2(i,j)=0.0d0
	dmut(i,j)=0.0d0
end forall
forall(i=1:n)
	ddelta(i,i)=1.0d0
	dphi0(2*n+1,i)=1.0d0
	dphi1(2*n+4,i)=c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4))
	dphi2(2*n+4,i)=-(c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4)))*(dsin(theta(2*n+5))**2.0d0)
	dphi2(2*n+5,i)=(c9-phi1(i))*2.0d0*dsin(theta(2*n+5))*dcos(theta(2*n+5))
end forall
forall(i=2:n)
	dc(n+i-1,i)=1.0d0
end forall

forall(i=1:m)
	dalpha0(i)=0.0d0
	dalpha1(i)=0.0d0
	dlambda(i)=0.0d0
	deta(i)=0.0d0
end forall
dalpha2=dalpha1
dalpha0(2*n)=1.0d0
dalpha1(2*n+2)=c9*2.0d0*dsin(theta(2*n+2))*dcos(theta(2*n+2))
dalpha2(2*n+2)=-dalpha1(2*n+2)*(dsin(theta(2*n+3))**2.0d0)
dalpha2(2*n+3)=(c9-alpha1)*2.0d0*dsin(theta(2*n+3))*dcos(theta(2*n+3))
deta(2*n+6)=1.0d0

forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:n)
	dmut(i,i)=1.0d0
end forall

gamma=phi0/(1.0d0-phi1-phi2)

gammat=diag(gamma)
igammat=diag(1.0d0/gamma)
forall (i=1:m,j=1:n)
	dgamma(i,j)=((1.0d0-phi1(j)-phi2(j))*dphi0(i,j)-phi0(j)*(-dphi1(i,j)-dphi2(i,j)))&
	           &/(1.0d0-phi1(j)-phi2(j))**2.0d0
end forall
dlambda=((1.0d0-alpha1-alpha2)*dalpha0-alpha0*(-dalpha1-dalpha2))/(1.0d0-alpha1-alpha2)**2.0d0

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)
sigmat=lamb(1)*cc+gammat
isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(1))))

jit=dot_product(epst(1,:),matmul(isigmat,epst(1,:)))

forall(i=1:n)
	sigeet(i)=(e_xi/(nu-2.0d0))*(dot_product(isigmat(i,:),epst(1,:))*&
	          & dot_product(isigmat(:,i),epst(1,:)))-isigmat(i,i)
end forall
siget=matmul(isigmat,epst(1,:))

dlnui=-(.5d0*n/(nu-2.0d0))+.5d0*jit*(e_xi/(nu-2.0d0)**2.0d0)&
&-.5d0*dlog(2.0d0)-.5d0*dpsi(.5d0/eta)+.5d0*e_logxi

grad=(e_xi/(nu-2.0d0))*matmul(dmut,siget)&
    &+.5d0*matmul(dgamma,sigeet)&
	&+.5d0*dlambda*((e_xi/(nu-2.0d0))*(dot_product(c,siget)**2.0)-dot_product(c,matmul(isigmat,c)))&
	&+lamb(1)*((e_xi/(nu-2.0d0))*dot_product(c,siget)*matmul(dc,siget)&
	&-matmul(matmul(dc,isigmat),c))&
	&+dlnui*(-1.0d0/eta**2.0d0)*deta



!t>1
do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0+alpha1*((fttm**2.0d0)+wttm)+alpha2*lamb(i-1)

	forall(ii=1:n)
	    metn(ii)=dot_product(igammat(ii,:),c)*dot_product(igammat(:,ii),epst(i-1,:))
	end forall
	matwtt=matmul(matmul(dgamma,ent),vec2((igammat.x.cc).x.igammat))

	dwttm=-(wttm*wttm)*(-(lamb(i-1)**(-2.0))*dlambda+2.0d0*matmul(dc,matmul(igammat,c))&
	      &-matwtt)
	
	dfttm=dwttm*dot_product(c,matmul(igammat,epst(i-1,:)))&
	     &+wttm*matmul(dc,matmul(igammat,epst(i-1,:)))&
		 &-wttm*matmul(dgamma,metn)-wttm*matmul(dmut,matmul(igammat,c))


	forall (ii=1:m,j=1:n)
		dgamma(ii,j)=dphi0(ii,j)+dphi1(ii,j)*(((epst(i-1,j)-c(j)*fttm)**2.0d0)&
		                                                                   &+(c(j)*c(j))*wttm)&
		           &+dphi2(ii,j)*gamma(j)&
				   &+phi1(j)*2.0d0*(epst(i-1,j)-c(j)*fttm)*&
				   &(-dmut(ii,j)-dc(ii,j)*fttm-c(j)*dfttm(ii))&
				   &+phi1(j)*2.0d0*c(j)*dc(ii,j)*wttm+phi1(j)*c(j)*c(j)*dwttm(ii)+phi2(j)*dgamma(ii,j)
	end forall


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)
	igammat=diag(1.0d0/gamma)


	sigmat=lamb(i)*cc+gammat
	isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(i))))

	jit=dot_product(epst(i,:),matmul(isigmat,epst(i,:)))

	forall(ii=1:n)
		sigeet(ii)=(e_xi/(nu-2.0d0))*(dot_product(isigmat(ii,:),epst(i,:))*&
		          &   dot_product(isigmat(:,ii),epst(i,:)))-isigmat(ii,ii)
	end forall
	siget=matmul(isigmat,epst(i,:))

	dlambda=dalpha0+dalpha1*((fttm*fttm)+wttm)+dalpha2*(lamb(i-1))&
	      &+alpha1*(2.0d0*fttm*dfttm+dwttm)+alpha2*dlambda


	dlnui=-(.5d0*n/(nu-2.0d0))+.5d0*jit*(e_xi/(nu-2.0d0)**2.0d0)&
	&-.5d0*dlog(2.0d0)-.5d0*dpsi(.5d0/eta)+.5d0*e_logxi

    grad=grad+(e_xi/(nu-2.0d0))*matmul(dmut,siget)&
		&+.5d0*matmul(dgamma,sigeet)&
		&+.5d0*dlambda*((e_xi/(nu-2.0d0))*(dot_product(c,siget)**2.0)-dot_product(c,matmul(isigmat,c)))&
		&+lamb(i)*((e_xi/(nu-2.0d0))*dot_product(c,siget)*matmul(dc,siget)&
		&-matmul(matmul(dc,isigmat),c))&
		&+dlnui*(-1.0d0/eta**2.0d0)*deta

end do

grad=-grad
end subroutine max_step_scr

end module studentscr